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We derive explicit spin and charge correlation functions of the N x N Hubbard 
model from a recently obtained weak-coupling analytic ground state The spin 

correlation function shows an antiferromagnetic behaviour with different signs for the 
two sublattices and its Fourier tranform is peaked at Q = (tt, tt). The charge correlation 
function presents two valleys at 45 degrees from the axes. Both functions behave in 
a smooth way with increasing A*'; the results agree well with the available numerical 
data. 



I. INTRODUCTION 



Let us consider the Hubbard model with hamiltonian 



H = Ho + W = tY^ cl^c,,^ + YuhrThri, U>0, 

cr {r,r') r 



(1) 



on a bipartite square lattice A = ^UBofA''xA'' sites with periodic boundary conditions and even A'^. 
Here a =t, J. is the spin and r, r' the spatial degrees of freedom of the creation and annihilation operators 
c' and c respectively. The sum on (r, r') is over the pairs of nearest neighbor sites and fira is the number 
operator on site r of spin a. The point symmetry is dv, the Group of a square besides, H is invariant under 
the commutative Group of Translations T and hence under the Space Group G = T ® C4„; (8) means the 
semidirect product. The presence of spin and pseudospin symmetries [1] leads to 5*04 Group [2] [3]; below, we 
shall work in the subspace of vanishing spin and pseudospin. We represent sites by r = {x, y) and wave vectors 



by k = {kx,ky) = -f{x,y), with x,y = Q 



,A'' — 1. In terms of the Fourier expanded fermion operators 
have Ho = X/k ^('^^''ko-^kCT with e(k) = 2t[coska: -t-cosfc^]. Then the one-body plane 
wave state cjj.^|0) = |k(j) is an eigenstate of Hq. 

The N X N Hubbard model at half filling is not elementary, even in the innocent-looking case of finite 
and small repulsion U. Indeed, weak-coupling expansions have long been known to be highly informative [4], 

/ 2N — 2 \ ^ 

[5] . The reason is that the trivial (7 = case is ( ^ ^1 times degenerate, and so even for relatively small 

lattices one has to solve a big secular problem to see how the interaction resolves the degeneracy in first-order. 
To deal with this problem, recently [6] [7] we have proposed a local formalism, based on diagonalizing the 
occupation number operators in the degenerate eigenstates of the kinetic energy Hq. We came out with an 
analytic singlet wave function |*I'[4^) which solves the secular problem and belongs to the ground eigenvalue 
of H . Under a lattice step translation it just flips spins {antiferromagnetic property). Further we proved that 



^C4v is the symmetry Group of a square. It is a finite Group of order 8 and it contains 4 one dimensional irreps, 
Ai, A2, Bi, B2, and 1 two-dimensional one called E. The table of characters is 
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it has vanishing momentum and its point symmetry is the same as the ground state symmetry established by 
Moreo and Dagotto [8]. 

We believe that 1^'^^;.) clearly deserves further study: although it is an eigenstate only at the first order in 
U , it must represent a good deal of the properties of the full ground state. Here we wish to show that the 

same local formalism that allows one to build l^*!^'^), is also suitable to bring out some physics. It is clear that 
the importance of analytic results is actually enhanced in the age of computers, since they can benchmark the 
numerical approximations. 

The correlation functions are a popular tool to understand and visualize the structure and the physical 
properties of a given many-body state. The definitions are: 

C?charge(r) = (*!;'^|nrno|*l;'^), (2) 

for the charge correlation function and 

G.pin(r) = (vE';,"^|Sr-S„|*;,%). (3) 

for the spin one. Here fir is the number operator and Sr is the spin vector operator at site r; the pedices 
denote the site at the origin. 

Studies [9] of correlation functions in the three-band Hubbard model have been aimed to the characterization 
of possible pairing mechanisms short after the discovery of high- Tc superconductivity [10]. Much of the early 
work dealt with the one-band model at strong coupling. Let us mention the exact diagonalization study by 
Kaxiras and Manousakis [11] on the vTo x \/T0 lattice, showing the antiferromagnetic order at half filling; 
the diagrammatic approach by Gebhard and VoUhardt [12] used the Gutzwiller ansatz, mainly for the Id 
chain. Correlation functions on larger lattices of the one-band [13], and, more recently, of the three-band 
Hubbard Model [14] [15] have been obtained by Quantum-Monte- Carlo methods. They have also been used 
to benchmark the self-consistent theory by Vilk et al. [16], which is basically a generalized Random-Phase 
Approximation. 

After summarizing the local formalism and the ground state solution in Section II, we go over in Section HI 
to a new picture by a particle-hole canonical transformation which is convenient to calculate the correlation 
functions. We derive the spin correlation function in Section IV and the charge one in Section V. The results 
are discussed and compared with available data in Section VI. 

II. THE GROUND STATE AT WEAK COUPLING 

In order to estabilish some notations we need to review the ground state formalism [6], [7]. Let Shf denote 
the set (or shell) of the k wave vectors such that £(k) = 0. At half filling {N^ particles) for ?7 = the Shf 
shell is half occupied, while all |k) orbitals such that e(k) < are filled. The k vectors of Shf lie on the square 
having vertices (±7r,0) and (0, ±7r); one readily realizes that the dimension of the set Shf, is \Shf \ = 2iV — 2. 
Since TV is even and H commutes with the total spin operators, 

= \Y.{nr^-nri), S+=^ct^Cri, = {S+)\ (4) 

r r 

at half filling every ground state of Hq is represented in the = subspace. Thus, Hq has 

degenerate unperturbed ground state configurations with = 0. 

It can be shown [6] that the structure of the first-order wave functions is gained by diagonalizing W in the 
truncated Htlbert space TL spanned by the states of N — 1 holes of each spin in Shf ■ In other terms, one solves 
a (2A'^ — 2)-particle problem in the truncated Hilbert space H and then, understanding the particles in the 
filled shells, obtains the first-order eigenfunctions of H in the full A'^^-particle problem. We emphasize that 
the matrix of Ho in Ti is null, since by construction 7i is contained in the kernel of Hq. 

The large set Shf breaks into small pieces if we take full advantage of the G symmetry. Any plane- wave 
state k belongs to a one-dimensional irrep of T; moreover, it also belongs to a star of k vectors connected by 
operations of C4v , and one member of the star has kx > ky > 0. We recall that any k £ 5^/ lies on a square 
with vertices on the axes at the Brillouin zone boundaries. Choosing an arbitrary k. £ Shf with kx > ky > 0, 
hence kx + ky = tt, the set of vectors Rik € Shf, where Ri € Cav, is a basis for an irrep of G. The high 
symmetry vectors k^ ~ (7r,0) and ks = (0, tt) are the basis of the only two-dimensional irrep of G, which 
exists for any N. If N/2 is even, one also finds the high symmetry wavevectors k = (±7r/2, ±7r/2) which mix 
among themselves under Cav operations and yield a four-dimensional irrep. In general, when k is not in a 



2 



special symmetry direction, the vectors iiik are all diflterent, so all the other irreps of G have dimension 8, 
the number of operations of the point Group dv 

Below, we shall need the number of these irreps. Since 8 times the number of eight-dimensional irreps + 
4 times that of four-dimensional ones -|- 2 for the only two-dimensional irrep must yield \Shf\ = '2N — 2, one 
finds that Shf contains A^e = — 2) irreps of dimension 8 if N/2 is even and No = 5('y ~ 1) irreps of 
dimension 8 if N/2 is odd. 

In this way, Shf is seen to be the union of disjoint bases of irreps of the Space Group. This break-up of 
Sh f enables us to define a real symmetry adapted one-body local basis which allows to carry on the analysis 
for any N. 

The one-body local basis is obtained by projecting onto the irreps of C4v the |k) states of Shf that belong 
to a given irrep of G. As already noted, kx = (vr, 0) and ks = (0,7r) belong to Shf and are the basis of a 
two-dimensional irrep of G. Let 



V2 



(|kA) + |kB)), liJB,) 



1 

71 



(|kA)-|kB» 



(5) 



be the first two real states of the local basis. As the notation implies, both axe simultaneously eigenvectors of 
the Dirac characters of C4v and carry symmetry labels; actually the symmetries are Ai and Bi because the 
two-dimensional irrep of G breaks into Ai ® Bi in Cav- In G these two functions merge into one irrep because 
the k states pick up phase factors from the translations. 

For even N/2, Shf also comprises the basis wave vectors ki = (7r/2,7r/2), k2 = (— 7r/2, 7r/2), ka = 
{tt/2, — 7r/2), k4 = (— 7r/2, —n/2) of the 4-dimensional irrep of G. This irrep breaks into Ai (B B2(B E in dv 
Letting I = 1, 2, 3, 4 for the irreps Ai, B2, Ex, Ey respectively, we can write down four more real local states 



(6) 



where O' is the 4x4 unitary matrix which performs the projections, namely, 



0' = 



(7) 



For N > 4, Shf also contains k vectors that are away from special symmetry directions. These form eight- 
dimensional irreps of G since i?ik are all different for all Ri £ C4,,. In other terms, any eight-dimensional irrep 
of G is the regular representation of Cav Thus, by the Burnside theorem, it breaks into Ai (B A2 (B Bi Q B2 ® 
E (B E, with the two-dimensional irrep occurring twice; these are the symmetry labels of the local orbitals we 
are looking for. Let k''"' = (fei'"', fel;"'') with fci™' > fcjj"' > be a wave vector of the m-th eight dimensional 
irrep of G and let 7?,, i = 1, . . . , 8 denote respectively the identity 1, the counterclockwise and clockwise 90 
degrees rotation , the 180 degrees rotation C2, the reflection with respect to the j/ = and a; = 

axis Ox, Oy and the reflection with respect to the x = y and x = —y diagonals cr^, a'y. We write down real 
local basis states as 



^o«|i^ikl'"l), 



(8) 



where O is the 8x8 unitary matrix 
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(9) 



Here, denoting by E' the second occourrence of the irrep 7 = 1, . . . , 8 is the A^, B2, E^, Ey, A2, Bi, E'^, 
E'y irrep respectively. 
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Now let us consider the following determinantal state 



m=l m=l 

for even and 

No No 

for odd N/2, with tr =|, J,. In Ref. [6] [7] we have shown that 

• |$af)ct is an eigenstate of W with vanishing eigenvalue {W — state). 

• Under a lattice step translation |$af)<t — > — |3>af)-<t- Therefore, it manifestly shows an antiferromag- 
netic order {antiferromagnetic property). 

r o] 

• Introducing the projection operator Ps on the spin S subspace, one finds that Ps\(i>AF)a = \^'Xf) ^ 
0, VS = 0, . . . , TV - 1. Then, {^af\W\^af) = Es^o {'^af\W\'^af) = 0, and this implies that there is 
at least one W = state of W in TC for each S. By the Lieb Theorem [1], only the singlet component 
I^af) belongs to the ground state multiplet of H at weak coupling (filled shells are understood, of 
course) . 

• has vanishing total momentum, is even under reflections, while the point symmetry is s or d for 
even or odd N/2, respectively. These are the correct quantum numbers of the interacting ground state 

at half filling [8] . 

• The ground state interaction energy per site is 

Eu = = J - U^;^. (12) 

Thus the linear term of the expansion of the energy per site in powers of U increases monotonically with 
N towards the infinite square lattice value U/4. 

Finally wc emphasize that in only 2N — 2 particles are antiferromagnetically correlated while in the 

strong coupling limit all the A'^^ particles show antiferromagnetic correlations. 

I^af) is an exact ground state of H for U ^ 0. Does this mean that it is the U ^ limit of the unique 
interacting ground state? We know that this is the case for the 4x4 and 6x6 square lattices, where we have 
numerical evidence that is the only singlet eigenstate in T( with vanishing eigenvalue. The correlation 

functions that we find below behave quite reasonably also for Af > 6 and this strongly suggests that 
continues to be a good approximation to the true ground state at weak coupling. A proof of the uniqueness 
of the vanishing eigenvalue of W in the singlet subspace of Ti. would be sufficient (although not necessary) 
to prove that. In the next Section we find further evidence to support this proposal: wc write I'l'^'p) in the 
particle-hole transformed picture and show that the corresponding Lieb matrix is positive semidefinite, as it 
should be for a genuine ground state [1]. 

III. THE GROUND STATE IN THE PARTICLE-HOLE TRANSFORMED PICTURE 

The unitary particle-hole transformation on a square N x N lattice and even TV reads 

r = {x,y). (13) 
This transformation maps the repulsive Hubbard model described in Eq.(l) onto the attractive one 



crT = (-r+^4T ' 



a {r,r') r 
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with fira' = dladrrr and -/Vct'*' = ^j.nr?. Letting {j'l'r)} be an orthonormal real basis of A'^^/2-particle states 
(that is, each l^fp) must be an homogeneous polynomial of degree N"^ /2 in the dl with real coefficients acting 
on the vacuum), we remind that the ground state at half filling 

|*[°1>= ^Lrir,|*riT>®l*r,i) (15) 

T1T2 

is such that the Lieb matrix irir2 is positive (or negative) semidefinite. 

In this Section we will perform the unitary particle-hole transformation in Eq.(13) on the ground state of 
Eqs.(10)(ll). We will show that the corresponding Licb matrix is indeed positive semidefinite; besides, it is 
already diagonal in the local basis. As a consequence, the ground state in the untransformed (c) picture is a 
pseudospin (as well as spin) singlet. 

Let ct, i = 1,... ,2N — 2, be the operators which create a particle in the i-th local state contained in 
\^AF)a, Eqs.(10)(ll). We write I^af)^ as 

I^Af),. =cl_^...c]v_i,^c]v,_^...ctjv-2,-a|0)- (16) 

It is clear that the first [last] N — 1 creation operators refer to the states of spin a [—a] in Eqs.(10)(ll). The 
singlet projection gives 

1 N-l 2JV-2 

|$ra ) = Ps=o\^AFh = E ^-^'sk E ^n-S,-,_, E ^i^-C} X 

k=N/2 zjv ,>..>ii=l ^ 3JV .>..>3i=iV ^ 

2 2 

Xc| ^ . . -cl^-i^fCN^l ■ ■ •C2jv-2,J.|0) + 

1 JV-1 2JV-2 

+;^{ E (-)'^"^^^ E E ^^-^J- 

k=N/2 iN_^>->il=i 3N._fc>-.>31=JV 

Xc| I • ■ • ^JV-l,i^]v,T ■ ■ ■ '^2iV-2,Tl'^) (■'■''') 

where = cj^Ci^i, = (S^)^, M is the normalization constant 

N/2 , .2 



and the p^'s are given by 



iV/2 + k-l 
M2 - 1 



N/2 
N/2-k 



Let 



1 ik r , 1 ik r , 



be the Fourier transformed operators of the site annihilation operators Cr and dr respectively. Prom Eq.(13) 
we get 

Cki = dki, CkT = dQ_kT> Q = (t^'"^)- (21) 

The ground state with the Fermi sea explicitly written is given by |^'J4J7) = I^af) ® |S) where |S) is the 
contribution from the filled shells: 

|E) = |ET)®|Ei), |E<,)= II cL|0). (22) 

e(k)<0 

Modulo an overall phase factor, the particle-hole transformation yields 

iSi)= n ^=^^110)= n 4iio)=iSi) (23) 

e(k)<0 e(k)<0 
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ist>= n '=kTio)= n dQ-^^iidU^)- (24) 

e(k)<0 e(k)<0 k 

Let rfj be the operator obtained substituing Ck with dk in the definition of c,. We note that e(k) < corresponds 
to e(Q — k) > 0. Then, the spin up filled shells state |E|) can be written as 

|St>= n 4t|0> =dI.T-4-i,T4,T-4jv-2,T n 4t|0> =4,t-4-i.t4,t-4jv-2,tI^T> (25) 

e(k)<0 e(k)<0 

and hence 

® |Ei) = dl^^..dj^_i_^<^^..4^_2_^|ET) ® |Ei). (26) 
The next step is to express Cia in terms of dia- By direct inspection one readly realizes that 

The above result implies that the raising operators S^' in the d picture are given by 

c+ ^ ( d^^d^i = D^ i^l,... ,N -1 . , 

' \-dndii = -Di i = N,...,2N-2 ■ 

These last three equations allow to rewrite the whole ground state |^'{^^) = I^aIt) ® |E) in the new picture: 

1 JV-l 2JV-2 



k=N/2 iN_^>..>ii=l ■'■N_fc>-->31=^ 



^'^i,T--'^'*f-i>T'^]v,i--4iv-2,j.4,t""4-i,t4,t ""4jv-2, 



1 JV-l 2JV-2 

k=N/2 «w:_fc>-->n=i ■?'^_fc>-->Ji=Ar 

|dAr,|..d2Ar-2,T'^i_|--iijv-i,t'^Ar,f'^2JV-2,Tl^T) ® l^^i) 
1 JV-l 2JV-2 

= 7^7^ E 5^ E E ^n-^L_.^-^--J X 



VAA — — — ' 

fc=JV/2 ijv >..>ii = ljjv ,>:>ji=N 
2 2 

x4,T-4jv-2,T4a-4jv-2,il^T) •X'lSi) + 

1 JV-l 2JV-2 

E E D,r-D.f.A-Dl^Jx 

fc=JV/2 i„_j^>..>ii = l j„_^>..>ii=Ar ^ 

x4,T--4-i,T4,i--4-i,il^T> ® (29) 
Therefore, the singlet ground state has the following form 

I*l4°l) = Y,wrVl^\t^)®Vl^\t^}, (30) 
r 

where F = {71, .., 7jv-i} with 1 < 71 < .. < 7jv-i < 2A'' — 2 and Z?^ = dJ^^ ..ci^jy_j . wr is the amplitude 
corresponding to the configuration F. If in F there are p indices between 1 and N — 1 and N — 1 — p indices 
between N and 2A'^ — 2, or viceversa, the amplitude wr is given by 

wr = -j=gN/2-p- (31) 

We conclude that the Lieb matrix is already diagonal in the local basis, and the nonvanishing diagonal 
elements arc TOr > 0. Thus it is positive semidefinite, which is consistent with its use as the ground state. We 
are now in position to calculate the correlation functions. 
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IV. THE SPIN CORRELATION FUNCTION 



In this Section wo will explicitly write down an exact analytic formula for the spin correlation function of 
the half filled Hubbard model in the limit of vanishing interaction. In particular we evaluate 

apin(r) = {v|/5i,jSr-So|*5i,), So = Sr=(0,0) (32) 

where Sr = {S^ , Sr,Sr) is the spin vector operator at site r with components 

S^ = l{S++S-), = ^^{S+ - S-), S^ = i(nrT-nri) (33) 



and 



Si- = cl^c,l, 4" = (5+)^ = ct^CrT. (34) 

Taking into account that is a singlet, one has (S^ So) = {Sr Sf^) = (SrSo) and hence 

Gspin(r) = ^ [{S-S+ + S+So)] (35) 

with = (Sj^Ljo 0) ^'^'^ (■ ■ ■) means the expectation value over the ground state l^'J^i?)- Since is a real 

linear combination of real basis vectors, G '"(r) = {S^ Sq) G 3i Vr and this implies 

G-+(r) = G-+(r)* = {SoS+) = (S+So) - 26rMSr). (36) 

Noting that is a translation invariant state, (Sq) = J2r(^^) = -^{S') = and hence 

G-+(r) = (S,+^,7)=G+-(r). (37) 

This last equation allow us to express the spin correlation function Gspin(r) in terms of G '"(r) only 

G.pi„(r) = ^G-+(r), (38) 

and the original problem is reduced to the calculation of G ^(r). In the following we will show that G ^(r) 
can be expressed in terms of three main contributions, two of which easily computable in the particle-hole 
transformed picture. Therefore, it is convenient to express G~"'"(r) in terms of the d's operators: 

G-+(r) = (cl^Cr^cl^coi) = {-r+^dl^dl^do^dol} = {-r+' y^r,wrA^^\Vr,idl^doiVljt^) x 

x{t^\Vr,A^doi-Dlji^) = 
= (-)'+" ^«'^lW^,G^l^,(r)^ (39) 

where, dropping the spin index, 

Grir,(r) = (ri|4rfo|r2), |r) = (40) 

Here and in the following coc = Cr=(o.())CT and doa = dr=(o,o)CT- Since the wr's are non-negative, Eq.(39) shows 
that G^^(r) is positive if r belongs to the sublattice A containing r = (0, 0) and negative otherwise. This was 
pointed out in Ref. [17]. All the information on the spin correlation function is enclosed into the site-dependent 
matrix elements Grir2(r). To evaluate them we write the annihilation operator dr as the sum of three pieces 

dr = ^Yl = ^'^W^d.ir) + ^ddr) + ^df(r) = P./di(r) + p[dj(r) + df(r)], 

JV ^ iV JV JV ^^^^ 



with 



di(r) = ^ e^'dk, ^^(r) = ^ — ^ e^'dk, df(r) = ^ — ^ e*"'dk 

V e(k)=0 V l*^! e(k)<0 V l"^! e(k)>0 (42) 
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and \Shf \ = 2iV - 2, |5| = i(iV^ - |5h/|). We observe that di(0) = di(r = (0,0)) belongs to Ai and that it 
can be written as a real linear combination of all the Ai-symmetric local annhilation operators of the local 
basis. By a unitary transformation on this Ai-subspace we may arrange that di (0) is the new di . Thus, from 
now on, the one-body local basis {d||0), i = 1, . . . , 2N — 2}, is such that the set of the yli-symmetric local 

states contains d\{0)\0) and 4|0) = 4(0)|0). 

Taking Eq.(41) into account one can express Grir2(r) as the sum of two terms: 

Grir,(r) = p^^/{ri|dl(r)di|r2) + {n\dl{r)d^\r2} = Pa/^Yr. W + p'^lflr, W, (43) 

with di = di(0) and d^ = df (0) = d^{r = (0,0)). Gp|p^(r) can be easily evaluated: 

4«P^(r) = (ri|4(r)ddr2)=5r,r, ^ ^ ^-^--^(El^dkJS) = 5r,r. ^ (44) 

e(ki)<0e(k2)<0 ' ' ' ' 

whore T(r) = J^^jj^j^,) e^* '' is the trace of the translation matrix in the Hilbert space spanned by the negative 
energies one-body plane-wave states. Taking into account Eqs.(43)(44), G~^(r) in Eq.(39) can be rewritten 
as 



(r) = (-) ^''{phf 2^ w;riWr2Gt,/^(r) +2p p^f-^ 2^WrG{,^\r) + p 2^Wr}. 



(45) 



By definition Wp = (^'[J'^l^'jJ^) = 1. To evaluate the diagonal matrix elements G{?/'(r) we need to use 

the antiferromagnetic property. In the local basis, the one-body translation matrix has an antidiagonal block 
form if a; -(- y is odd and hence a diagonal block form otherwise. Therefore di(r) can be expanded as 

{I]^^^ii(i')di, ii(r) e K x + y even 
(46) 
Yh=n^ ti{r)di, ti{r)€^ x + yodd 

and G[?/' (r) becomes 

Gi!'/l(r) = (rldl(r)d.lr) = I ^^W^^- ^ + ^-7 (47) 

where 71 is the first index of the configuration F (we remind that F — {71, . . . ,7Ar-i} with 1 < 71 < . . . < 
7JV-1 < 2N — 2). Substituting this result into Eq.(45) we see we need to evaluate ^p WpJi-y^. This can be 
done by observing that I^JJi?) can be rewritten in the c picture as 

|$ra) = |$ra)_|$[oi) (48) 



with 



JV-2 JV-1 2JV-2 

Vat 



k=0 ik>..>ii=2jk>-->jl=N 



JV-1 2JV-2 



|<') = 47E(-)V. E E ^t-^s^;;-^;;<i---4-M4,T---4jv-2.Tio) (so) 



k=0 ik>..>ii=2ji,>..>ji=N 



and 



f _ j 9N/2-k k = 0,... ,N/2-l 

- \ gk+i-N/2 k = N/2,... ,N-1 ■ ^^'^ 

All the configurations contained in ]$|°') are such that in the particle- hole transformed picture ^i-^j = 1. On 
the other hand, all the configurations contained in are such that in the particle-hole transformed picture 

(5i7i = 0. Therefore 
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The second term in Eq.(45) is totally determined onco we know ti{r). By definition 



ti(r) = (0|dUr)di|0) = ^1^, 



(53) 



where Th/(r) = 



(k) = 



is the trace of the translation matrix in the Hilbort space spanned by the e(k) = 
one-body plane- wave states. In the last equalility of Eq.(53) we have used Eq.(42). By noting that 7h/(r) 
vanishes any time x + y is odd, one obtains for G '"(r) the following result 



G-+(r) = i-r+^ptf wr.wr.G^^flirf + T(r)[T(r) + T,,(r)]} 



(54) 



where we have used Eqs.(47)(52)(53). 

In order to make this result more explicit, we perform the sum in the first term of Eq.(54). It can be easily 
calculated coming back to the original c picture. Indeed 



(55) 

where is defined in Eq.(50) and T'(r) is the translation operator by r: T'^(r)ciT'(r) = ci(r). As usual a 

is given by the same expression which defines di, with dk — » Ck. Therefore, according to the new expression 
of di — di(0), we have Cj — ^ y^,.f|,)-n '^k and this is why ci — ci(0) = ci(r = (0,0)) shows up in 
Eq.(55). Finally wo observe that the spin-dependent filled Fermi sea jEo-) can contribute only a phase factor 
corresponding to its momentum; since jSo-) has vanishing momentum the phase factor is exactly 1. 
The explicit evaluation of X{r) is deferred to Appendix A. Here we report the final result 



X(r) = 



(-)" 



X + y oven 



A + B X (4iV - 4) x + y odd 



(56) 



where A and B are two A''-dependent constants. Eventually, substituing this last result in Eq.(54) and taking 
into account Eq.(38), we get the full analytic expression of the spin correlation function 



T(r)\T(r)+T,f(r)] + i ^ + B^T^f{r) x + y even 
^ [r) + lhf[r)\ + S ^ + g ^ (4iV - 4) x + y odd 



(57) 



In this form it is not hard to show that independent of the numerical value of the two constants A and B the 
sum rule 



^G.pi„(r) = 



(58) 



holds. Indeed, let us consider the identities 

^(-l)^+«T(r)2 = ^ ^ (_^)X+yg-i(k+k').r ^ ^^_i(k+k'+Q).r. ^gg) 

r r e(k),e(k')<0 £{k),e(k')<0 r 

since the r summation yields A''^ times a S function, and the 5 function is never satisfied (if e(k) < then 
e(Q — k) > 0), one finds that 



^(-l)-+«T(r)2 = 0, 

r 

r 

e*" = ^("^k.o + 5k,Q), 

where A is the sublattice with sites having x + y even, one gets 



and, similarly 



On the other hand, since 



(60) 



(61) 



(62) 



(63) 



and hence Eq.(58). 
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V. THE CHARGE CORRELATION FUNCTION 



The charge and spin correlation functions are closely related. Let nlr*^' = hri ± "hri', so, fii^' is the number 
operator, while hi ' is twice the z component of the spin on the site r. Then, the charge correlation function 
is given by 



Gcha.ge(r) = (*51.|nWnW|*J,°l,) , 
while the spin correlation function can be written as 



where h, 



[±1 = ^[±1 



nlJ((,o). Let 



|*2U = |*[°1>-|*M), 



(64) 



(65) 



(66) 



with l^'IJ'') = (S) be a suitable decomposition of the singlet ground state wave function, see Eq.(48). 
Exploiting the invariance of (^'CT'|n[*'nQ*'|^'|^j) for simultaneous flips of a and a' we get 



^charge ('^) — 2 



U+lUT,[01 



l^i[+luT/[ol^ 



(67) 



J^spin(r) — 2 



(<>|nHr(r)nH|<>) - (<>|nHf (r)nH|*M) 



where T(r) is the operator of the translation by r, such that 



(68) 



(69) 



and T{r)\^^^p) = \^^af) ^as been used. The action of h\^^ on the state j^'i"') can be easily evaluated. We 
can express Cr as the sum of three operators like in Eq.(41) 

cr = PhfCi{v) + p[cj(r) + Cf (r)] (70) 

where ci, and Cf are defined as di, and d| in Eq.(42), but dk must be substituted with Ck. Then we get 



«-o*'l*r) = (noT ±n-oi)|*|'^) = {pis + ± p^) + pphf{c\ ci^^ ± c\^^ca) + P^dc^:-^ ± dc^i 



(71) 



4*'!*™} = (^OT ± ft0i)|*f') = (P^ ± Plf ± P^) + PPh}[c\ .^C^^ ± 4^Cl,|) + P^(4^CJT ± 4^ca 



1*^ 



(72) 



Hence (not ^ ^oi)|*<T) can be expressed as a linear combination of five orthogonal states. By means of these 
two last equations one gets 



(*^">^*'T(r)n™|*^"') = [{pis + p'± py + 2p^(S|4^Cf,r(r)4^cjT|S>]($^"'|T(r)|$^"') + 

ct 



p'pis{{^^°^\cl,f{vM^f){nc-,,f^^^^^ (73) 



and 



(*[°>|nWT(r)n™|*^"'> = [±(p^, + p^ ± p^)^ + 2p^(E|4,Cf,T(r)4^cjT|E>]($i"' |T(r)|$^"') 



U±1|,T,[01 



[0]| 



(74) 



since {Y,\(^^^c^^T{r)c^^^c^a)\Y,) does not depend on cr. The number of scalar products can be further reduced: 
we recall that |3'J^i^) = l^f) — l^j"') is an eigenstate of the total momentum with vanishing eigenvalue: 



l = ($51.|T(r)|cE.[,«^)=2 



<I.M|f(r)|$^"')-(<I.i"JlT(r)||<I.^"J) 



and hence 



($f=|f(r)|$^°i) = «i|T(r)|$i"i)-l 



[0]\ 



(75) 



(76) 
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Substituting Eq.(76) into Eq.(74) and then subtracting Eq.(73) term by term yields 



, 2 2 
+P Phf 



+ [{i±{-mpi,+p'±pr] (77) 



+ 



and so, since Phf + 2p^ = 1, 



C?charga(r) = 1 + |apin(r) + ptf[l - 4y(r)], 



(78) 



where 

y(r)s($l°l|f(r)|$[°l). (79) 
We postpone to Appendix B the exphcit calculation of ^(r). Here we limit to present the final results 



y(r) 



D + ExThf{r)'^ x + yeven 



~ + IsJ^]^ ^ ] D + Ex{4N-4) x + yodd 



(80) 



As for the spin correlation function one can easy verify that independent of the numerical value of the two 
iV-constants D and E, the sum rule 



y^<3charge(r) = 



(81) 



is satisfied. 



VI. RESULTS AND DISCUSSION 



Most of the available data on the half filled Hubbard Model on a square lattice refer to the 4x4 cluster, 
see for example Ref. [5]. In the left hand side of Figure 1 we report a classical representation of the spin 
correlations in the 4x4 lattice: the length of the lines is proportional to the absolute value of the correlation 
function, and the sign is positive for the lines going up. This representation was adopted in Ref. [5] and 
our result is identical to that reported there, which was obtained by second-order perturbation theory on the 
computer. 




FIG. 1. Left. The spin correlation function between the origin (empty circle) and the other sites; the length of the 
lines is proportional to the absolute value of the correlation function, and the sign is positive for the lines going up. The 

on-site value is reduced by a factor 0.4 for graphical convenience. Right. The spin correlation function in real spa<;e for 
the 4x4 model, along a counterclocliwise path from the origin (empty circle, see the inset). 
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More data [18] [19] on the 4x4 cluster were obtained by Fano, Ortolani and Parola by exact diagonalisation 
augmented by an intensive use of Group theory techniques. In the right hand side of Figure 1 we report the 
spin correlation function in real space, Gspin(r) in Eq.(57), along a triangular path. Although our results are 
exact for U 0, remarkably the trend is quite the same as the one reported in Ref. [18] for (7 = 4. An overall 
factor of 4 depends from the definition of the spin operators in [18], lacking the usual 1/2 factor. 

The analytic expression of the spin correlation function in Eq.(57) agree with the important Shen-Qiu-Tian 
[17] theorem, which has been extended to finite temperatures quite recently [20]. This theorem states that 
the spin correlation function must be positive on one sublattice and negative on the other. This applies to 
the results for the 4x4, 6x6, 8x8 and 10 x 10 clusters as well, as shown in Figure 2. Essentially, the 
Shen-Qiu-Tian property is a consequence of the positive semidefinite ground state Lieb matrix, and we have 
explicitly verified this property in Section III above. 




6 1 




The Fourier transform of the spin correlation function 

api„(k) = J2 e* ""apin(r) (82) 

r 

is shown in Figure 3; the ticks on the x axis correspond to the points F = (0,0), P = (7r,0), Q = (tt, tt) and 
r again, in fc-space. The trend is seen to converge rather quickly to a characteristic shape which is strongly 
peaked in the Q direction. 
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rPQr rPQr 




rpQr rpQr 

FIG. 3. The Fourier transform of the spin correlation functions in clusters of varoius sizes. The ticks of the x axis 
correspond to the F, P, Q and F, as usual. 

The charge correlation function in real space shows characteristic structures with two intersecting channels 
at 45 degrees from the axes as exemplified in Figure 4 for the 10x10 case; at the intersection the correlation 
function presents a narrow hole. Similar trends are observed for the other clusters, although the intensity of 
the corrugation declines with increasing cluster size. 




FIG. 4. Real space representation of the charge correlation function in the 4x4, 6x6, 8x8 and 10 X 10 clusters. 

The Fourier transformed charge correlation function is dominated by a delta function at F = (0, 0) resulting 
from the almost constant distribution in real space. In Figure 5 we have removed that delta. The figure 
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represents Gchargo(k) = e* '"Gchargc(i') along the path F, P, Q and F (see figure caption for details) for 
the 10 X 10 square lattice. Already at this cluster size Gchargc (k) shows a very similar trend to its asymptotic 
(A'' oo) shape. 
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10.4 
o 



oL| \ \ \i 

r p Q r 

FIG. 5. The Fourier transform of the charge correlation function in the 10 xlO cluster; a delta function at the F 
point is removed. The ticks of the x axis correspond to the F, P, Q and F, as usual. 



In conclusion, we have obtained explicit analytic expressions for the spin and charge correlation functions 
using a weak-coupling ground state wave function |*5f) of the half filled Hubbard Model on a square lattice. 
We compared our analytic results with the numerical data available in the litterature. They always agree well; 
remarkably, provided that U < t/N^, our predictions are good approximations to the exact diagonalization 
results. 

As far as the non half-filled system is concerned, the same local formalism can be used to calculate the 
correlation functions of the doped Hubbard antiferromagnet, but first-order perturbation theory is not enough 
to single out a unique ground state. However, in the half filled ground state there are 2N — 2 particles in the 
e = shell that do not have double occupation; therefore, doping the system with two holes, one obtains a 
first-order ground state provided that a pair is annihilated belonging to the shell Shf, this must be a VF = 
pair. Since there are W = pairs belonging to different irreps of the Space Group, the many-body ground state 
which is formed by annihilating the pair also has components of different symmetries. For each symmetry, we 
shall have different correlation functions to compute. Actually, we need second-order perturbation theory to 
resolve the degeneracy. This was done in Ref. [21] in the special case A'' = 4. On the other hand, the problem 
becomes trivial when the shell at the Fermi surface is totally filled since the non-interacting ground state is 
unique. 

As pointed out in Ref. [22], one can use standard perturbation theory to calculate the correlation functions 
order by order in U; in the thermodynamic limit they can be expanded as an asymptotic series. 



APPENDIX A: EVALUATION OF X 



To evaluate X{r) — {^^^^\c\ ^ci^^T{r)c\ .^ci.||$'j°') we need to consider separately the cases of even and odd 
X + y. In the first case, let T{r) be the block-diagonal translation matrix in the local basis: 

f EL7r(rK.ct,jO), i = l,...,A-l 
r(r)4|0) = J2 r(r).,,4l0) = <^ . (Al) 

[ E^^^'nr)..^4i0), i = A,...,2A-2 

We have 



2JV-2 JV-1 2JV-2 



r(r)ci_^ci,j$r') J2 En ^w-- n 

ai..ajv-l = l /31'-/3n-1=-'V a = l 6=JV 
N-2 JV-1 JV-1 

E(-)'^ E E ^X,--5i,^iS,-^..5,-.^cl^,^..ct^_^,^4,,..ct^_^,j0) (A2) 



X 

fe=0 ifc>..>ii=2it.>..>ji=l 



and hence 
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iV-1 2JV-2 N-1 2N-2 P 



k X 



ai..c«jv-i = l ;3i ..;3]v-i=iV a=l 



/c=0 



Af-1 JV-1 



^{ E E (0|c~-i.i-ci,i5;;,..5-,5r5X^..5X^5+ct^,^..ct^_^,jO>}x 

2iV-2 N-1 

E E <0|c2iv-2.T..c;v.T54..5+5^.^..5^,^ct^,^..ct^_^,^|0)}. (A3) 

nfc>..>ni=JV jj.>..>ji=l 



Taking into account that the annihilation operators in the second and third row of the above equation are all 
different (in particular their indices are 1, . . . , iV — 1 in the second row and N, . . . , 2N — 2 in the third row), 
a great simplification takes place: 



<fc>..>il=2 



<ji>-->il=l,7^ai 



(A4) 



(A5) 



J)s>->J1 = 1 



3k>-->jl=N 



where e is the totally antisymmetric tensor with N — 1 indices, while e/3i../3^_i = e^j_jv+i..;ajv-i-Jv+i and the 
dots mean that we are neglecting other terms whose contribution to the scalar product is zero. On the right 

hand side of Eq.(A4) the summation indices ik > ■■ > ii run in the interval {1, ... , AT — 1} in such a way that 
none of them is equal to ai. Using Eq.(A4), the second row of Eq.(A3) yields 



m/j>..>mi=2 ij(.>..>ii=l,^ai 



-ai..aiv-i 



while using Eq.(A5), the third row of Eq.(A3) yields 

2N-2 2N-2 

£f}i..0N-i E E (Olc2iv-2,T--Cjv,T'§4--'^"i'^3-^--'^A'=Af,T--4jv-2,T|0) 



nfc>..>ni=JV jfc>..>ji=JV 



These two last results allow to rewrite X{r) as 



— e/3i../3jv_i 



N-1 



^« = P^+-5Ci,i(r)T(r)i,i 



where, using the convention j ' , ) = for a binomial coefficient with negative down entry. 



A — '^hfT ■r^JV-2 r2 I N — 1 
^ ^ 7} Z-,k=0 J k 



N -3 
k-1 



(A6) 



(A7) 



(A8) 



(A9) 



1 ■^N-2 ,2 t N - 

77 2^k=o \ 



N -2 
k 



N -3 
k- 1 



(AlO) 



and Ci,i(r) is the (1,1) algebraic complement of the matrix T(r) whose determinant is equal to one for even 
X + y. The (1,1) algebraic complement can be expressed in terms of the (1,1) element of the matrix T(r): 
Ci,i(r) = Tt(r)i,iDet[T(r)] = T(r)i,i (since T{r)ij € K). Next one has to recognize that T(r)i,i is, by 
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definition, equal to ti(r) (see Eq.(46)) whose analytic expression is given in Eq.(53). Therefore, any time x + y 
is even we have 



(All) 



On the other hand, for odd x + y the translation matrix in the local basis is antiblock diagonal, that is 

E^=^'nrkc4|0), i = l,...,N-l 

E^="i' r(r)i,^ct |0), i = N,...,2N-2 



(A12) 



7=1 



and hence 



JV-l JV-l 



ai..aN-l=JV/3i../3jv-l=l a=l 



N-1 N-1 



x^{-)''fk E E '^^ifc--'^ii'^«i'^/3:(fc--'^/3ji'=li,r-^L-i,i4i,T--4iv-i,Ti°)- 



(A13) 



ifc>..>il=2 jfc>..>ji=l 



Let us consider the fe-th term of this state. One can check by direct inspection that the only non vanishing 
contribution in the scalar product with the state c\ ^ci,||$j°') comes from the {N — k — 2)-th term of the 
corresponding expansion, see Eq.(50). Hence 



JV-l JV-l 



w = -j^ E EH ^('■)'''"<' n ^('■^'''''^ E f^^''-"- 



2 X 



i E 

mN-k-2>--~ 

2JV-: 

i E 



ai..aiv_i=JV/3i../3jv-l = l o=l 
JV-l 



mfj_k_2>-.>mi=2jj,>..>ji = l 

2JV-2 JV-l 

> E (0|c2^-2.T-CJv,T5+„_,_,..5isX^..SX^S+cl^,^..cl„_^,jO)}. (A14) 



njv-Jc-2>-->"l='V <ji>-->ii=2 



Analogously to the case x + y even, the particular structure of allows the following simplification 



jfc>-->ii=i 



jfc>..>ji=i 



(A15) 



<ic>--><l=2 <ji>..><l=JV,ySai 



(A16) 



where eai..aj^_i = eQi-jv+i...ajv-i-JV+i and the sum on the right hand side of Eq.(A16) means that the 
indices ik > ■■ > ii run in the interval {A^, . . . , 2N — 2} in such a way that no of them is equal to ai. All the 
neglected terms contribute nothing to the scalar product. 
Using Eq.(A15) the second row of Eq.(A14) yields 



e/3i../3N_i E E (0|cjv-ia--cia'5'm„_fc_2--5,„i<5i ^■^..S'j^cl .f..c]v_i,||0) = 

' N-2 

— e/3i../3jv-l 



while using Eq.(A16) the third row of Eq.(A14) yields 

2JV-2 2JV-2 

Sai-aw-l E E (0|C2JV-2,T--Cjv,t5'+^_^_2..S+^S,^..S,+ SJ^c]v,^..4jV-2,J0) 

njV-fe-2>-->"l=^ i)c> ->il=JV,#ai 



(A17) 



-ai..ajv-l 



N-2 
k 



(A18) 
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Substituing these results in the expression for X(r) we get 



where we have taken into account that Det[T(r)] = — 1 for odd x + y and the constant C is given by 



^ \Sh.fl^ f f ( 



iV-2 
k 



(A19) 



(A20) 



In the last equality of Eq.(A19) we have used C = A + B x (AN — 4) which is a direct consequence of the sum 
rule for the spin correlation function, see Eq.(58). 



APPENDIX B: EVALUATION OF Y 

Here we show that Y{v) = {<I>^"'|T(r)j<l?'j°') has the form shown in Eq.(80) and we derive the explicit values 
for the two constants D and E. As for X{r) we will first consider the case x + y even and thereafter the case 
X + y odd. Making use of Eq.(Al), we get 



2N-2 N-l 



b,l3b X 



ai..aiV_i=l/3i../3iV-l=JV a=l 



E(-i)^^ E E 



b=N 



t + 



1 l^3k" "I'T" "N-l.T 0N,l" 02N-2 



.JO) 



(Bl) 



ik>--->ii=2jk>--->jl=l 



The fc-th term in the sum of Eq.(Bl) gives non- vanishing scalar product only with the fc-th term in Eq.(49) 

and hence 



2N-2 N-l 



=jf n En ^w-- n E /' 



k X 



b=N 



k=0 



N-l 



ai..o:]si-i=l Pi..I3n-1=N a = l 
N-l 

x{ E E (0|cJv-i,T--Ci.T.?m,--5+iS-^..S-^c+^_.r..ct^_^^^|0>} X 

7nif.> ..>m,i—2 ii^> . .>ii=2 
2N-2 N-l 



(B2) 



Now we use the fact that the indices ai,. . . , un-i must be all different and within the range {!,... , A'' — 1} 
otherwise the scalar product vanishes. This means that 



il^>..>il=2 



N-l 

'^ii--'^ifeCl,f ^]v-1,t|0) + ■ ■ ■ 

ifc>..>n = i,7^ai (B3) 



where the neglected terms do not contribute to the scalar product. In the second row of Eq.(B2), ci,| commutes 
with all the raising spin operators whatever are the values of the A; indices rnt > .. > mi in the range specified 
by the sum. This implies that ii cannot be 1. On the other hand cl,^ commutes with all the lowering spin 
operators and hence no one of the indices rrik > . . . > mi can be ai otherwise the corresponding term vanishes. 
Hence the term in the secon row gives 



E E Sai--aN- 

mfe>..>m.i=2,7^ai ii^>..>ii_=2,^ai 



^{0\CN-l,r-Cl,tSmk-SmiSi^..Si^c\^^..cl,_j^^^\0) 



= [ Yl (B4) 

ij, >...><! =2, T^ai 

If = iV — 2 the sum is zero except for ai = 1 since the indices ik > ■ ■ ■ > ii don't have space to run. Hence 
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= l5fc,JV-2<5l,ai + (1 — 5k,N-2) 



ifc>..>il=2,7^ai 



N-3 



= (1 - 5fe,Ar-2) ( 



Sk,N-2 + (1 — Sk,N-2){ 



N-2 



k 

N-3 



(B5) 



while the third row of Eq.(B2) yields 

2Ar-2 2JV-2 



n^>..>ni=N jk>-.>ji=N 



e/3i../3jv_i 



iV- 1 



(B6) 



as can be verified by using the total antisymmetry of each homogeneous polynomial in the raising spin oper- 
ators. 

Therefore for even x + y one can write 



Y(r) = \ ^ - i?Ci.i(r)r(r)i,i = i + ^^[D + Ex T,j (rf 



where D and E are two iV-dependent constants given by: 



D 



(B7) 



(B8) 



'=^|(7V-l)/^_ 



J.2 / N -1 
2 + l^k=l Jk 



N-2 
k 



N-3 
k 



For odd x + y we make use of Eq.(A12). Then, the action of T'(r) over |3>|.°') gives 

2Ar-2 Af-1 N-1 2N~2 

f (r)i3.M) = ^ ^ En ^^W"."'' n ^('■)''-''^ ^ 



ai..aiv_i=Ar/3i../3N-l=l o=l 



b=N 



fc=0 <fc>...>ii=2 jfe>...>3i=l 



In the scalar product with the j^i"') state, only the terms with the same number of up (or down) spins in 
the first N — 1 (and hence in the last — 1) local states will survive. Let us consider for example the fe-th 
term of the sum in the second row of Eq.(BlO); it contains states where k of the last N — 1 local states have 
spin down and k of the first N — 1 local states have spin up. This term has non-vanishing scalar product only 
with the (A'^ — k — l)-th term of the sum in the definition of Eq.(49). In particular this implies that the 

terms where the first and the last N — 1 local states have all the spins aligned do not contribute to the scalar 
product. Hence 

2iV-2 JV-1 JV-1 2Ar-2 JV-2 



(B9) 



(BIO) 



ai ..c<_v_i=JV /3i../3iv_i=l a=l 

N-1 



'"JV-fc-l>-->'ni=2 jfc>-.>il=l 
2JV-2 JV-1 
X-l > ^ (0|C2Ar-2,i-CjV,i'S', 



«{ E \U|C2Ar-2,i-CjV,i-3n„_fc_ 

nisi-k-i>..>ni_=N i^>..>il=2 

Let us consider the term in the second row. For a given choice of /3i../?jv-i one finds 

JV-1 JV-1 

E ^4-^44ii-4«-iiio) = ^'3i-/3«-i E ^it-5;,cia..4-i,iio) + . 
jfc>->ji=i jfc>->ji=i 



(Bll) 



(B12) 
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where the missed terms do not contribute to the scalar product. Since the Ci,| annihilation operator coninmtes 
with all the raising spin operators coming from the sum over niN-k-i > .. > mi, ji is constrained to be 1 
for all non vanishing contributions. Still for jk > .. > j2 fixed there is only one choice for rriN-k-i > ■■ > mi 
to have non- vanishing result. In particular, the possible results for a given choice of jk > ■■ > 32 and 
ruN-k-i > ■■ > mi are or 1. Hence the term in the first square bracket yields 



JV-l 



N-l 



mjv _ Jt _ 1 > ■ ■ >iTii =2 jfe > . . > J2 =2 

f 

-Sqi..qN-l y k_l 

A similar trick can be used for the term in the third row of Eq.(Bll). We get 

N-l 2N-2 



ik ••'-^n''iV,T ••''2JV-2 



J0) + .. 



(B13) 



(B14) 



ifc>..>il=2 



ife>..>il=JV,7^ai 



Since the c|^^ ^ creation operator commutes with all the lowering spin operators coming from the sum over 
ik > ■■ > ii, one of the indices njv-fc-i > .. > ni is constrained to be ai and the third row of Eq.(Bll) can 
be rewritten as 



2N-2 



2N-2 



£ai..aiv-i ^ ^ 

nfi_k-l>-->n2=N,^ai if.> . .>ii=N 



(0|C2Ar-2,i--CjV,i>S'nN-fe- 



il^N,T-^2N-2 

N-2 



,tIo> = 

(B15) 



Substituting these results in Eq.(Bll) one obtains 



Yir) = ^ = i - ^[D + Ex{4N-A)] 



(B16) 



where we have taken into account that Det[T(r)] = — 1 for odd x + y and the constant Fn is given by 



— 77 — l^k=l JkjN-k-l 



N-2 
k-l 



N-2 
k 



(B17) 



In the last equality of Eq.(B16) we have used Fn ~ \Shf\'^ /4: — D — E x {4N — 4) which is a direct consequence 
of the sum rule for the charge correlation function, see Eq.(81). 
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